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In this study we provide the first numerical demonstration of the effects of turbulence on the mean Lorentz force and 
the resulting formation of large-scale magnetic structures. Using three-dimensional direct numerical simulations (DNS) 
of forced turbulence we show that an imposed mean magnetic field leads to a decrease of the turbulent hydromagnetic 
pressure and tension. This phenomenon is quantified by determining the relevant functions that relate the sum of the tur- 
bulent Reynolds and Maxwell stresses with the Maxwell stress of the mean magnetic field. Using such a parameterization, 
we show by means of two-dimensional and three-dimensional mean-field numerical modelling that an isentropic density 
stratified layer becomes unstable in the presence of a uniform imposed magnetic field. This large-scale instability results in 
the formation of loop-like magnetic structures which are concentrated at the top of the stratified layer. In three dimensions 
these structures resemble the appearance of bipolar magnetic regions in the Sun. The results of DNS and mean-field nu- 
merical modelling are in good agreement with theoretical predictions. We discuss our model in the context of a distributed 
solar dynamo where active regions and sunspots might be rather shallow phenomena. 
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1 Introduction 

Turbulence effects generally refer to the occurrence of corre- 
lations between velocity, temperature, and/or magnetic fields 
at small scales. A typical example is turbulent viscosity, 
which results from the spatial exchange of turbulent eddies 
characterized by velocity correlations. This leads to the dis- 
sipation of energy at small scales. 

However, there is also the possibility of additional (e.g. 
non-diffusive) turbulence effects, as is perhaps best known 
in mean-field electrodynamics and dynamo theory. Here one 
models the effects of the mean electromotive force, i.e. the 
turbulence effects of velocity and magnetic field fluctua- 
tions, on the evolution of the mean field. This can lead to 
the occurrence of the famous a effect, in addition to turbu- 
lent magnetic diffusion, turbulent diamagnetic velocity, etc. 
(Moffatt 1978; Krause & Radler 1980). Another example 
is the A effect in rotating anisotropic hydrodynamic turbu- 
lence, which can lead to the occurrence of differential rota- 
tion in cosmic bodies such as the Sun (Riidiger 1980, 1989; 
Riidiger & Hollerbach 2004). In that case the relevant cor- 
relations come from the mean Reynolds stress tensor and its 
dependence on the local angular velocity. 

A related example is the combined Reynolds and Max- 
well turbulent stress tensor and its dependence on the mean 
magnetic field. The first analytic calculations of the depen- 
dence of the turbulent Reynolds stress on the mean mag- 
netic field in the framework of the first-order smoothing ap- 
proximation were performed by Radler (1974) and Riidiger 
(1976). Later, also the combined effects of the turbulent 
Reynolds and Maxwell stress tensors were considered (Klee- 




orin et al. 1989, 1990; Riidiger & Kichatinov 1990). It was 
noticed that this can lead to a local reduction of the total tur- 
bulent pressure and hence to the possibility of self-induced 
concentrations of large-scale magnetic fields (Kleeorin et 
al. 1989, 1990, 1996; Kleeorin & Rogachevskii 1994; Ro- 
gachevskii & Kleeorin 2007). Such a process may play an 
important role in the formation of sunspots and active re- 
gions in the Sun. It may be complementary to the mag- 
netically induced suppression of the turbulent energy flux, 
which would lead to further cooling and hence a further con- 
centration of the structures (Kitchatinov & Mazur 2000). 

The current leading explanation for the formation of sun- 
spots is related to the emergence of deeply rooted magnetic 
flux tubes (Parker 1955, 1982, 1984). Such flux tubes are 
generally believed to be produced and 'stored' near the bot- 
tom of the convection zone (Spiegel & Weiss 1980). The 
storage of magnetic fields and the formation of flux tubes in 
the overshoot layer near the bottom of the solar convective 
zone was investigated in a number of publications (see, e.g., 
Spruit 1981; Spruit & van Ballegooijen 1982; Schiissler et 
al. 1994; Moreno-Insertis et al. 1996; Tobias et al. 2001; 
Tobias & Hughes 2004). However, in order that the tubes 
retain their basic east-west orientation during their ascent 
over many pressure scale heights, the magnetic field must be 
strong enough (Choudhuri & D'Silva 1990) and is estimated 
to be around 10 5 G at the bottom of the convection zone 
(D'Silva & Choudhuri 1993). Such fields would be up to a 
hundred times stronger than the equipartition value, which 
is one of several other arguments that have led to the idea 
that flux emergence of dynamo-generated fields might in- 
stead be a shallow phenomenon (Brandenburg 2005; Schat- 
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ten 2009). Such a scenario appears also compatible with so- 
lar subsurface flows, as inferred from local helioseismology 
(Zhao et al. 2001; Kosovichev 2002). In particular, Zhao et 
al. (2004) and Hindman et al. (2009) find the presence of 
converging flows around active regions at radii as large as 
100-200 Mm. It appears that these convergent flows might 
actually be the source of the formation of active regions and 
perhaps sunspots rather than a consequence (Parker 1979a; 
Hurlburt & Rucklidge 2000). Of course, in the immediate 
proximity of individual spots one observes outgoing flows. 
Those are probably superficial, less than 1-2 Mm deep, and 
appear to be due to the dynamical effects of magnetocon- 
vection in an inclined magnetic field of the penumbra (e.g., 
Thomas et al. 2002; Heinemann et al. 2007; Rempel et al. 
2009; Kitiashvili et al. 2009). 

The goal of this paper is to investigate the effects of 
turbulence on the mean Lorentz force by means of direct 
numerical simulations (DNS) for forced turbulence and to 
study the instability of a uniform large-scale magnetic field 
in an adiabatically stratified layer by means of mean-field 
numerical modelling based on parameterizations both of an- 
alytic formulae by Rogachevskii & Kleeorin (2007) and the 
results of DNS for forced turbulence. In order to study the 
essence of the effect, we make several simplifications by ne- 
glecting the energy equation, i.e. the specific entropy is as- 
sumed to be strictly constant in space and time. In the mean- 
field numerical modelling we neglect the suppression of tur- 
bulent magnetic diffusivity and turbulent viscosity, and omit 
correlations with density fluctuations. Nevertheless, the mean 
density is allowed to evolve fully self-consistently accord- 
ing to the usual continuity equation. 

2 Turbulence effects on mean Lorentz force 

Throughout this paper we adopt units for the magnetic field 
where the vacuum permeability is equal to unity, i.e. the 
magnetic pressure is given by \B 2 . 

2.1 General considerations 

We use the equations of mean-field magnetohydrodynam- 
ics (MHD). These equations are obtained by averaging the 
original MHD equations over small-scale fluctuations. This 
technique is best known in the case of the induction equa- 
tion (Moffatt 1978; Krause & Riidler 1980). In this study we 
are mainly interested in effects of turbulence on the mean 
Lorentz force. Let us consider the momentum equation, 

M pU = -^ (1) 

where 

Ily = p UiUj + S l3 ( P + \B 2 ) - BiBj - aij (2) 
is the momentum stress tensor, U and B are the velocity 
and magnetic fields, p and p are the fluid pressure and den- 
sity, Sij is the unit Kronecker tensor, <7y = 2pvSij is the 
viscous stress tensor, with 

(3) 



being the traceless rate of the strain tensor, and v is the 
kinematic viscosity. Ignoring the turbulent correlations with 
density fluctuations for low-Mach number turbulence, the 
averaged momentum equation is 

9 ~ (4) 



where 11^ = ELy + is the mean momentum stress ten- 
sor split into contributions resulting entirely from the mean 
field (indicated by superscript m) and those of the fluctuat- 
ing field (indicated by superscript f). The tensor ELy has the 
same form as Eq. (O, but all quantities have now attained 
an overbar, i.e. 

tl% = - P u^ 3 + (p + \w) - b~b 3 - a l3 . (5) 

We emphasize here that p is just the mean gas pressure and 
cfy is the average of the microscopic viscous stress tensor, 
<7y. The contribution from the fluctuating fields, in turn, 
is split into contributions that are independent of the mean 
fields (and hence isotropic and proportional to Sij) and con- 
tributions which depend on the mean fields, 

Tt Pto*« <• (6) 

Here, pto 1S the turbulent pressure in the absence of a mean 
magnetic field and <j°f = cry + afj quantifies the turbu- 
lent viscosity, ofj = 2~pv t Sij, and the additional effects of 
the mean magnetic field on the effective stress tensor. Both 
terms, pto$ij — <yjf, result from the fluctuations u = U — U 
and b = B — B of velocity and magnetic fields, and are 
given by 



PtrA 



= pUiUj + \5rjb 2 ~ bibj. 



(7) 



J V u ij 

In the absence of a mean magnetic field, the turbulent back 
ground pressure is 



Pto = iK + jP u o ( 8 ) 
where the subscripts on b 2 and u 2 indicate values in the 
absence of the mean magnetic field. Magnetic fluctuations 
b are generated both by small-scale dynamo action and by 
tangling of the mean magnetic field by velocity fluctuations. 
On the other hand, the velocity fluctuations also have two 
contributions, those that depend on the mean magnetic field 
and those that do not. 

Following Rogachevskii & Kleeorin (2007), the part of 
the effective stress tensor that depends on the mean mag- 
netic field is parameterized as 



o%(B) = —q s BiBj + \5 l3 q p B 2 



(9) 



W^ + U, hl )-\8 l3 V-U 



where q s and q p are functions of the mean field, B, and 
the magnetic Reynolds number, Rm. Eq. (O implies that 
the effective mean Lorentz force that takes into account the 
effects of turbulence, can be written as: 

pT M = -|V[(1 - q P )B 2 } + B ■ V [(1 - q s )B] . (10) 

The detailed analytic expressions for q s (B) and q p {B) have 
been given by Rogachevskii & Kleeorin (2007). The asymp- 
totic formulae for the nonlinear functions, q p (B) and q s (B), 
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are given below. For this purpose we define (3 = B /B cq , 
where B cq = (pit 2 ,) 1 / 2 is the equipartition field strength. 

For very weak mean magnetic fields, 4/3 <C Rm -1 ' 4 , 
<7 p and q s are approximately constant and given by 

«p(/3) 4 



45^ 
9.C9) = ^ (1 



■ 91nRm) 
41nRm) 



[l + 5|ln(4/3)| 
In(4/3)| 



32 /3 2 ] 
f 12/3 2 



for Rm~ 1/4 < 4/3 < 1 we have 
16 
25 

32 fiiT,r4^i + — 

30 

and for strong fields, 4/3 ^> 1, we have 

q p {(3) = 1/6/3 2 , a.09) = ^/48/3 3 . 

(Rogachevskii & Kleeorin 2007). 

In Sect. 12. 21 we present DNS evidence that the functions 
q p (B) and q s (B) are positive, indicating the possibility of 
a reduction of the effective Lorentz force, i.e., a decrease 
of the effective magnetic pressure and magnetic tension in 
small-scale turbulence. 

2.2 DNS of turbulence effects on mean Lorentz force 

In order to study turbulence effects on the mean Lorentz 
force and to determine the functions q p (B) and q s (B) from 
DNS, we consider forced turbulence in a periodic three- 
dimensional domain in the presence of an imposed uniform 
magnetic field, say Bq = (Bo, 0, 0). We determine q s and 
<7 p from Eqs. (Q and (O for i = j = x, 



)B 2 



0) 



(ID 



Pto + \{b 2 x -b 2 y - bl) - pu 2 x = [\q p - q s 
and i = j = y, 

p w + ^-W x -M)-pl4 = ±q p Bl (12) 

where p t o is given by Eq. {8). First, we determine p t o from a 
simulation with Bq = 0. Then, we use Eq. dTZb to determine 
q p (B). Finally, to determine q s (B) we subtract Eq. (fTTl i 
from Eq. (TT2l . i.e. 



2^ - 



(13) 



In order to determine separately the effects of the mean field 
on the turbulent Maxwell and Reynolds stresses we also 
compute their respective contributions q^ + q^ = q p and 
9^ + 9? = Qs- In the current case where there is no small- 
scale dynamo action we have 

qf = m-W-¥ z )/B 2 , <g = 2(p t0 -p^ y )/B 2 , (14) 

In Fig. Q] we show the results from the simulations of 
forced turbulence in the presence of an imposed magnetic 
field. Here we have also plotted the total turbulent energy, 
Et = Ek + Em where 
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Fig. 1 Bq dependence of the normalized turbulent energy 
Et — Et /Eq, where Eq is the value of Et for Bq = 0, 
together with the contributions from kinetic and magnetic 
energies, Ek = Ek/Eq and Em = Em/Eq, (upper panel) 
as well as the coefficients q p , q£, q^ 1 , and q s obtained from 
DNS for Re = 180 and R m = 45 (second and third panels). 



are the energy densities of velocity and magnetic fluctua- 
tions, respectively. The simulations were performed with 
the Pencil Cod^I which uses sixth-order explicit finite 
differences in space and third order accurate time stepping 
method. They are similar to the DNS of Haugen & Bran- 
denburg (2004). The forcing function consists of plane non- 
polarized waves with an average wavenumber fcf = 1.5 fci, 
where ki — 2ir/L is the smallest wavenumber that fits 
into a 3D periodic domain of size L. We use an isothermal 
equation of state with constant sound speed c s . The forc- 
ing strength is arranged such that the turbulent rms velocity, 
u rms , is below c s . In all our runs the maximum turbulent 
Mach number, u rms /c s , is around 0.2, so that compressibil- 
ity effects are still weak. 

The fluid Reynolds number based on the forcing wave 
number k{ = 2?r//f, is Re = u lms /vkf w 180 and the 
magnetic Prandtl number is v/rj = 0.25, so the magnetic 
Reynolds number, R m = u Tms /i~ik{ « 45, is just small 



E K =pu 2 /2, E M = b 2 /2 



(15) 



http://www.nordita.org/software/pencil-code/ 
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enough so that no small-scale dynamo is excited. [Follow- 
ing Haugen et al. (2004) and Schekochihin et al. (2005, 
2007), the critical value of R m is between 70 and 80 for 
this value of the magnetic Prandtl number]. Note that the 
magnetic Reynolds number, Rm, used in Sect. 2.1 is based 
on forcing scale If = 2n/kf, and thus Rm = 2nR m . In 
all cases we normalize the field strength in terms of the 

1 12 

equipartition value, £> oq = p u rms . The number of mesh 
points is 64 3 . 

As follows from our DNS study, both functions q p (B) 
and q s (B) are positive and exceed unity for weak fields. The 
error bars are obtained by calculating averages over each of 
three equally long intervals of the full time series and tak- 
ing the largest deviation from the full averages. For very 
small values of Bq / B cq , the turbulent fluctuations dominate 
over the effects from shredding the imposed field, which in- 
creases the error bars. Nevertheless, the results are in agree- 
ment with the prediction of Rogachevskii & Kleeorin (2007). 
Note also that both theory and simulations suggest that q p > 
2q s . 

2.3 Physics of turbulence effects on Lorentz force 

The physics of the effect of turbulence on the large-scale 
Lorentz force is as follows. The equation of state for isotropic 
turbulence is given by 
1 2 

Pto = ^E M + -E K , (16) 

[see, e.g., Landau & Lifshitz 1975, 1984, and Eq. ©], where 
Pto is the total (hydrodynamic plus magnetic) turbulent pres- 
sure, and Ek and Em are defined by Eq. $15[ . The total en- 
ergy density Et = Ek + Em of homogeneous turbulence 
with a mean magnetic field B is determined by the equation 

^=/T-^+m(VxB) 2 , (17) 
at t 

where To is the correlation time of the turbulent velocity 
field in the maximum scale If of turbulent motions, It is 
the energy source of turbulence, rj t is the turbulent mag- 
netic diffusion. For a given time-independent source of tur- 
bulence It and for t ^> to the total energy density of the 
turbulence reaches a steady state 

Et = const = tq It , (18) 

where we neglect a small magnetic source 7ft (V x B) 2 of 
the turbulence [that is of the order of 0{l 2 /H B ), and Hb 
is the characteristic scale of the mean magnetic field spatial 
variations]. The approximate constancy of Et is compati- 
ble with DNS, where we found only a small decrease (20%) 
for strong (equipartition strength) mean fields (see the up- 
per panel of Fig. [1). The reason for the departure is possibly 
a dependence of To on B. However, this decrease only en- 
hances the modification of the Lorentz force by turbulence. 

Equations ( fT6b and ([T8l imply that the change of turbu- 
lent pressure Spto is proportional to the change of the mag- 
netic energy density 5Em, in particular Sp t o — —(1/3)SEm 
(because SE^ = —SEm)- Therefore, the total turbulent 



pressure is reduced when magnetic fluctuations are gener- 
ated. 

For a non-zero large-scale mean magnetic field B, the 
change of the magnetic energy density SEm, is proportional 
to B 2 . Therefore, the total turbulent pressure is given by 
Pt = Pto — \% B 2 , where p t o is the turbulent pressure 
in a flow with a zero mean magnetic field. The coefficient 
<7 P > when magnetic fluctuations are generated. Now we 
examine the part, P c s(B), of the total pressure, P to t = 
p + P t + \B 2 , that depends on the mean magnetic field, B, 
i.e., we consider the effective mean magnetic pressure that 
takes into account the contribution of turbulence, P c ff (B) = 
|(1 — <7 p ) B 2 . We study the case when p 3> B 2 /2, and 
therefore, the total pressure P to t is always positive, while 
the effective mean magnetic pressure P e g(B) may be neg- 
ative when q p > 1, 

The modification of the mean Lorentz force can result 
in the excitation of a long wavelength instability even in an 
initially uniform mean magnetic field in a density stratified 
layer (Kleeorin et al. 1990, 1996; Rogachevskii & Kleeorin 
2007). Indeed, the growth rate of the instability for the per- 
turbations perpendicular to both the gravity field g and the 
mean magnetic field B is given by 



A = |h/(1-?P)(^-1 



(19) 



where c a — B/p^ 2 is the Alfven speed, H p 1 = Vlnpo|> 
Hg 1 = |VlnP|, and we neglected for simplicity the dis- 
sipation processes due to turbulent viscosity and turbulent 
magnetic diffusion. For q p > 1 the instability is excited 
when Hp < Hb, i.e., it occurs even in an initially uniform 
mean magnetic field. 

The mechanism of the instability can be understood as 
follows. An isolated magnetic tube moving upward becomes 
lighter than the surrounding plasma since the decrease of 
the magnetic field in it, due to expansion of the tube, is 
accompanied with an increase of the magnetic pressure in- 
side the tube. This increase, due to the fact that the effec- 
tive magnetic pressure is negative, leads to a decrease of 
the density inside the tube and to the appearance of a buoy- 
ancy force. It results in the further upward displacement of 
the flux tube, i.e. it causes the excitation of the instabil- 
ity. The instability causes the formation of inhomogeneous 
magnetic structures. The energy for this instability is sup- 
plied by the small-scale turbulence. In contrast, the free en- 
ergy in Parker's magnetic buoyancy instability, that is ex- 
cited when Hp > Hb, is drawn from the gravitational field 
(Newcomb 1961; Parker 1966). The growth rate of Parker's 
magnetic buoyancy instability is determined by Eq. ( fT9b for 
9p = 0. 

Magnetic buoyancy in astrophysics applies usually to 
two different situations (see, e.g., Priest 1982). The first cor- 
responds to a problem described by Parker (1966, 1979b) 
and Gilman (1970a, 1970b). They considered the instability 
of a stratified continuous magnetic field and did not invoke 
a magnetic flux tube. The other situation was considered 
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by Parker (1955), Spruit (1981), Spruit & van Ballegooijen 
(1982), Ferriz-Mas & Schussler (1993), and Schussler et al. 
(1994), who studied buoyancy of horizontal magnetic flux 
tubes. 

In the present study we investigate the first situation, i.e., 
we study the large-scale instability of a continuous (diffu- 
sive) magnetic field in small-scale turbulence. This instabil- 
ity is caused by turbulent velocity and magnetic fluctuations 
and leads to the formation of large-scale magnetic structures 
(see Sect. 3). 

3 Mean-field numerical modelling 

In order to understand in more detail the appearance of mag- 
netic structures from this large-scale instability we consider 
now numerical solutions of the mean-field MHD equations 
in a density stratified layer. We apply a new mean-field model 
which includes the effect of turbulence on mean Lorentz 
force. We solve the evolution equations for mean velocity 
U, mean density p, and mean vector potential A in the form 



f)lj 

— = —U ■ VU — Cg V lnp 

dt s ' 



M 



'K.tot, 



d lnp 

dt 



U Vlnp- V-U, 



QA _ 

— = UxB-( Vt + rf)J , 



(20) 



(21) 



(22) 



where ^F M is given by Eq. ( [Tol l, and ^K.tot 
with p-Fk = V • cr K and with p-FVisc — V • cr, so that 

^K,tot - (vt + v) (V 2 C7 + VV • U + 2SV lnp) (23) 

is the total (turbulent and microscopic) viscous force and S 
is given by Eq. <[3j. 

We consider two- and three-dimensional models of an 
isentropically stratified atmosphere, where the gravitational 
potential is written as &(z) = (z — z^g, so the gravity 
vector g = — V<I> = (0, 0, —g) is constant. We arrange the 
initial profiles of density and sound speed such that they 
take given reference values at z = 0, i.e. p = Po and c s = 
c s o at z — 0. This implies that z,^ = (3/2)c% /g. Our initial 
profiles therefore obey 



P/Po = (c s /c s0 ) 3 



3 



(24) 



The local density scale height is H p = / g, and the density 
scale height at z = is H p o — c^ Q /g. In Fig. [2] we show the 
vertical dependence of the initial density. 

We allow for the presence of an imposed field in the y 
direction, Bq = (0, Bo, 0). The total field is then written as 

B = Bo + VxI, (25) 

so the departure from the imposed field is expressed in terms 
of the mean magnetic vector potential A. 

On the upper and lower boundaries we adopt stress-free 
boundary conditions for velocity, i.e. U x>z = U V:Z = U z = 
0, and a normal-field condition for the magnetic field, i.e. 



10 



-10 



-6 -4 



-2 



Fig. 2 Density stratification in a model with z top = H p q. 
The dashed lines indication the reference value with p = po 
at z = 0. 



A z = for the vector potential. Here, commas denote par- 
tial differentiation. No boundary condition for the density 
is required. Again, all computations have been carried out 
with the Pencil Code. 

In setting up our model we define the wavenumber fcf of 
the energy-carrying eddies. This relates the turbulent mag- 
netic diffusivity to the rms velocity via rj t — u rms /3fcf . This 

means that the ratio of our non-dimensional field strength to 

1/2 

turbulent diffusivity, i.e. £>o/c s op to -qt/c^Hp, is given 
by 3k{Ho times Bo/B eq . 

In this paper we approximate q p and q s by simple profile 
functions, 



/ 2 W 

Qp = 9 P o I 1 - - arctan 



2 W 
q s = q s o I l ~ ~ arctan ^ 



(26) 



(27) 



Following Rogachevskii & Kleeorin (2007), we also define 
<5 p = 1 — q p and Q s = 1 — q s . Correspondingly, we de- 



fine the coefficients Q 



p0 



1 - q p0 and Q s0 = 1 



<?s0- 



B , 



By — B = 0, corresponding to A x ^ 



A 



For our fiducial model we use QpQ — 2C2so — 20, corre- 
sponding to q p0 = 21 and q s0 = 11. The resulting magnetic 
pressure is shown in Fig.[3]for their analytic theory and the 
result from DNS is shown in Fig.Hl together with the corre- 
sponding fits. For our fiducial model we choose furthermore 
— — 1/2 

B p = B s = 0.1 c s oPg • For the imposed field strength we 

1/2 

choose -Bo/ c soPo = 0.01 and for the turbulent magnetic 
diffusivity we take rj t /c s oH p = 0.01. As discussed above, 
this means that Bo/B cq = 1/3 if we assume k[Ho = 1. 

3.1 Magnetic structures in two-dimensions 

In this paper we consider both two-dimensional and three- 
dimensional solutions. We begin with two-dimensional mod- 
els in the xz plane with an imposed field in the normal (y) 
direction, B = (0,B ,0). 
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B/B eq 

Fig. 3 The effective mean magnetic pressure P m (B) — 
(1 — q p )B /B p determined by Rogachevskii & Kleeorin 
(2007) - solid line, and by the model described by Eq. (f26t 

— 1 /2 

- dashed line (B p = 0.21 c s oP and q p0 = 4). 




B/B eq 

Fig. 4 Same as Fig. [3] but from simulation (dotted 
line). The solid line shows a fit [Eq. d26j >l with B p = 

0.022 CsoPq (corresponding to B p /B cq — 0.18) and 

q P o = 21. 

In Table Q] we present a summary of some exploratory 
runs where we list the nondimensional growth rate A = 
\Hp/j]t, as well as the nondimensional saturation values 
of rms mean velocity and mean magnetic field, U lms = 
U lms H p0 /r]t and B rms = B Ims H p0 /r) t pQ . These experi- 
ments show that decreasing Q p o and Q s o lowers the growth 
rate and the saturation values, while increasing the degree 
of stratification (e.g., increasing zt op = zt op /H p o from 1 to 
1.2, corresponding to an increase of bottom to top density 
ratio from 108 to 237) enhances them. Likewise, increasing 
the strength of the imposed field enhances the growth rate 
and the saturation values, while increasing the magnetic dif- 
fusivity lowers them. These results are in agreement with 
Eq. 03. 

In Fig. [5] we compare the evolution of the rms values of 
velocity and magnetic field for two different stratifications. 
As discussed above, increasing the amount of stratification 



Table 1 Summary of non-dimensional run parameters to- 
gether with the resulting non-dimensional growth rates A as 
well as the non-dimensional saturation values of rms veloc- 
ity and magnetic field. The tildes indicate non-dimensional 
quantities, as explained in Sect. 13.11 Our fiducial run is 
Run C. 

-Bo fjt Ztop QpO QsO A f/sat -Bsat 
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0.01 


0.01 
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0.000 
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0.01 
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0.013 
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0.01 


0.01 
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-20 
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0.013 
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0.01 
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-20 


-10 
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0.025 
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0.01 
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0.001 


0.001 
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0.013 
0.026 
0.010 
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0.0001 
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Fig. 5 Growth of the rms value of mean velocity and 
mean magnetic field for two runs with different degree of 
stratification. 



1.000 




_0.010 

B a /c aa pl /z 



0.100 



Fig. 6 Growth rate as a function of Bq, keeping all other 
parameters as for the fiducial Run C. 



increases the growth rate. The scaling of the growth rate 
with the strength of the imposed field is shown in Fig. [6] 

Next, we consider the structure of the resulting magnetic 
field. In Fig.|7]we show "meridional" (xz) cross-sections of 
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Fig. 7 Early evolution of magnetic field in the y direction 
(color coded) together with velocity vectors in the xz plane. 
Note the spontaneous production of flux structures. 

the magnetic field at three different times during the early 
phase where the magnetic field just begins to saturate. The 
horizontal wavelength of the pattern is about 10 H p . As time 
goes on, the structures of enhanced magnetic field inten- 
sify. The decrease of the effective pressure makes them even 
heavier which explains their subsequent descent. 

At later times new structures can form near the surface. 
In Fig. [8] we show an example during the fully saturated 
phase of the instability, where one sees the emergence of a 
new patch that is then swept to the larger one and merges. 

3.2 Magnetic structures in three-dimensions 

In order to demonstrate the three-dimensional nature of the 
instability we extend the domain in the y direction, so both 
x and y are between ±10H p o, and —10 < z/H p q < 1. 
These simulations are otherwise similar to those in the two- 
dimensional cases. In Fig. [9] we show visualizations of the 
field at three characteristic times. Note in particular that 
the wavelength of the pattern in the y direction (parallel to 
the field) is 3-4 times shorter than that in the x direction 
(perpendicular to the field). Again, the instability begins to 
emerge first near the surface where it develops magnetic 
structures which begin to sink downwards. Viewed from 






B/B 



Fig. 8 Later evolution of magnetic field for the same run 
as in Fig. [7] Note the mutual merging of flux structures. 



above, one sees the emergence of what looks like multiple 
bipolar regions. 

In Fig. [10] we show a horizontal cross-section from an- 
other simulation, where the horizontal extend is only half as 
much as before. This figure is suggestive of the formation of 
bipolar structures. The plane has been rotated by 90° such 
that the y direction points now from left to right. The black 
horizontal bar gives the density scale height at the depth of 
the cross-section, which is about twice the value H p o at the 
reference depth. 

4 Conclusions 

In this study we have demonstrated in DNS the effects of 
turbulence on the mean Lorentz force. This effect is quanti- 
fied by determining the relevant functions q p (B) and q s (B) 
that relate the sum of the turbulent Reynolds and Maxwell 
stresses with the Maxwell stress of the mean magnetic field. 
Using three-dimensional simulations of forced hydromag- 
netic turbulence with an imposed field, we confirm that the 
function q p (B) is positive and can reach values much larger 
than unity for B/B cq <C 1. This thereby reverses the sign 
of the effective magnetic pressure P e g(B) associated with 
the mean magnetic field, which then becomes P e g(B) = 
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A. Brandenburg et al.: Large-scale magnetic flux concentrations 






rection, so the y direction, corresponding to the azimuthal 
direction when applied to the Sun, points from left to right. 
BjB The black horizontal bar gives the density scale height at the 

depth of the cross-section, which is about twice the value 
H n n at the reference depth. 

0.0115 p * 



0.0000 
-0,0115 
-0.0229 

O 

Fig. 9 Visualizations of the magnetic field at the early 
saturation phase (t = 700 H p o/c s o = 7 H^/rjt left), at an 
intermediate time (t = 1000 H p0 /c sQ = 10 Hjp/rft left), 
and a later time (t = 1300 H p0 /c s0 = 13 !!%/■)% left). 
Note the broad similarity of field in the xz plane with the 
two-dimensional cases. Note that the wavelength of the pat- 
tern is shorter in the y direction than in the x direction. In 
the final time much of the magnetic field structures have 
sunk beneath the surface, leaving only a few isolated bipo- 
lar structures at the surface. 

i(l - q p )B . We find that the function q s (B) that deter- 
mines the modification of magnetic tension, also reaches 
values much larger than unity, but its value is less than half 
the value of q p and the error bars are larger. This work has 
also demonstrated explicitly the possibility of a large-scale 
instability of the full system of mean-field MHD equations. 
Finally, our study has revealed the presence of spatial struc- 
tures arising from the instability that might be associated 
with the formation of bipolar magnetic regions and perhaps 
also sunspots when applied to the Sun. 



The effects of turbulence on the large-scale Lorentz force 
may also be important in applications to the solar torsional 
oscillations by changing the mean magnetic tension B ■ VB 
(Rudiger et al. 1986; Rudiger & Kichatinov 1990). More 
specifically, these turbulence effects may be critical for ex- 
plaining the narrow structure of the observed solar torsional 
oscillations (Kleeorin et al. 1996). 

Clearly, there are several possibilities of improvement 
that might make the model more realistic and eventually 
suitable for confrontation against observations. Most impor- 
tant is perhaps the fact that our reference values B p q and 
B s o in the quenching profiles (1261 and (l27l i. as well as the 
coefficients g p o and q S Q in these profile functions, are cur- 
rently kept constant. It will be more realistic to make them 
vary with depth, because density and turbulent rms veloc- 
ity vary with depth. Another extension of the model would 
be to allow for the possibility that the magnetic field to be 
generated by a mean-field dynamo rather than relying on an 
imposed field. 

On the more theoretical side, there is a need to further 
verification of the essential physics of the negative mag- 
netic pressure effect. In particular, one must wonder why 
the effects of this instability have not yet seen in DNS. A 
likely possibility is the lack of sufficient scale separation. 
Only now realistic simulations are beginning to be large 
enough to encompass sufficiently many cells in all three di- 
rections. The importance of sufficient horizontal extent was 
already emphasized by Tao et al. (1998), who find clear ev- 
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idence of a segregation into strongly magnetized and weak 
magnetized regions, a phenomenon that might be closely 
related to that reported here. A particularly promising ap- 
proach might be to generalize the direct simulations dis- 
cussed in the present paper to the case with vertical den- 
sity stratification such that the setup becomes similar to the 
mean field models that we also studied in this paper. The 
turbulence would then still be driven by a forcing function. 
Another possibility is to study this effect with turbulent con- 
vection instead of forced turbulence. This is particularly in- 
teresting, because theoretical predictions of Rogachevskii & 
Kleeorin (2007) suggest that the modification of the effec- 
tive Lorentz force will be even stronger in turbulent convec- 
tion. 
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